*Code replicates and produces main figures in the manuscript.

*************************************************************
*Figure 1: Hydraulic Fracturing Potential and Production - Rystad Prospectivity Index (RPI)

*Figure 1b requires restricted data from Enverus
*************************************************************
*Panel A
use "$tmp/tquartilemap.dta",clear
qui count if tquartile1 == 1 
loc n1 = r(N)
qui count if tquartile1 == 0
loc n0 = r(N)
spmap map_category using "$tmp/cntyfips2000coords.dta" ///
, id(id) ///
title("") ///
clmethod(unique) fcolor(white gs10 gs6) ///
ndf(white) ndo(white) /// 
legend(position(4) size(*1.3) order(2 "No Fracking Potential" 3 "Shale Play Counties (`n0')" 4 "Top-Quartile Counties (`n1')"))
graph export "$out\fig1A.pdf", as(pdf) name("Graph") replace

*Panel B
use "$tmp\fracking_production.dta",clear
twoway (line allprod_hwell0 year, ///
lcolor(blue) lpattern(dash) lwidth(medium)) ///
(line allprod_hwell1 year, ///
lcolor(blue) lpattern(solid) lwidth(medium)) ///
(bar n_plays year if !missing(n_plays), ///
barwidth(0.8) color(gray%30) yaxis(2)) ///
(scatter n_plays year if !missing(n_plays), ///
msymbol(none) mlabel(bar_label) mlabposition(12) ///
mlabcolor(black) mlabsize(small) yaxis(2)), ///
xlabel(1990(5)2020) ///
ylabel(0(500)2500, axis(1)) ///
ylabel(0(1)4, axis(2) nolabels noticks) ///
yscale(axis(1) range(0 2500)) ///
yscale(axis(2) range(0 4)) ///
ytitle("Oil and Natural Gas Extracted: Millions of BOE", axis(1)) ///
ytitle("", axis(2)) ///
xtitle("") ///
legend(order(1 "Other Shale Play Wells" 2 "Top-Quartile Wells") ///
position(6) cols(2) region(lcolor(none))) ///
graphregion(color(white)) plotregion(lcolor(black))
graph export "$out\fig1B.pdf", as(pdf) name("Graph") replace

*************************************************************
*Figure 2: Earnings and Employment Effects by Gender
*************************************************************
use "$tmp/final_playonly.dta", clear
foreach var in earns emptotal {
	forval i = 1/2 {
		qui reghdfe ln`var'_gen`i'_a1499 ddtquartile1 $rhs $main
		lincom ddtquartile1
		loc beta1: di %4.3fc r(estimate)
		loc beta = trim("`beta1'")
		loc se1: di %4.3fc r(se)
		loc se = trim("`se1'")
		loc p = r(p)

		reghdfe ln`var'_gen`i'_a1499 alpha_* $rhs  $main 
		eventstudy gen`i' 0 "`beta'" "`se'" `p' o -.03 .03 .06
		// Four panels: A=earns male, B=earns female, C=emptotal male, D=emptotal female
		if "`var'" == "earns" & `i' == 1 {
			local panel "A"
		}
		if "`var'" == "earns" & `i' == 2 {
			local panel "B"
		}
		if "`var'" == "emptotal" & `i' == 1 {
			local panel "C"
		}
		if "`var'" == "emptotal" & `i' == 2 {
			local panel "D"
		}
		graph export "$out/fig2`panel'.pdf", replace 
	}
}



*************************************************************
*Figure 3: Overall Mortality by Gender
**Requires Restricted CDC Mortality Data
**Variable naming in code (See readme for more details)
**cdrdc0_gen0_eth00_a2564, crude overall death rate (dc0), both genders (gen0)
**cdrdc0_gen1_eth00_a2564, crude overall death rate (dc0), male (gen1)
**cdrdc0_gen2_eth00_a2564, crude overall death rate (dc0), female (gen2)
*************************************************************
use "$tmp/final_playonly.dta", clear
forval i = 0/2 {
	qui reghdfe cdrdc0_gen`i'_eth00_a2564 ddtquartile1 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)

	reghdfe cdrdc0_gen`i'_eth00_a2564 alpha_* $rhs  $main 
	eventstudy gen`i' 0 `beta' `se' `p' 
	// Add panel letters: A=overall, B=male, C=female  
	local panel = cond(`i'==0, "A", cond(`i'==1, "B", "C"))
	graph export "$out/fig3`panel'.pdf", replace 
}

*************************************************************
*Figure 4: Internal vs. External Causes of Death
**Requires Restricted CDC Mortality Data
**Variable naming in code: 
**cdrdc30_gen0_eth00_a2564, crude internal death rate (dc30), both genders (gen0)
**cdrdc31_gen1_eth00_a2564, crude external death rate (dc31), both genders (gen0)
*************************************************************
use "$tmp/final_playonly.dta", clear
*internal causes of death
foreach d in 30 31 {
	qui reghdfe cdrdc`d'_gen0_eth00_a2564 ddtquartile1 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)
	reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* $rhs  $main 
	if `d' == 30 {
		loc panel "A"
		loc ymin = -40
		loc yint = 20
		loc ymax = 40		
	}
	else if `d' == 31 {
		loc panel "B"
		loc ymin = -20
		loc yint = 10
		loc ymax = 20				
	}
	eventstudy dc`d' 0 `beta' `se' `p' o `ymin' `yint' `ymax'
	graph export "$out/fig4`panel'.pdf", replace 
}

*************************************************************
*Figure 5: Internal Causes: Cardiovascular vs Other Internal Mortality
**Requires Restricted CDC Mortality Data
**Variable naming in code: 
**cdrdc4_gen0_eth00_a2564, crude cardiovascular death rate (dc4), both genders (gen0)
**cdrdc32_gen1_eth00_a2564, crude non-cardiovascular internal death rate (dc32), both genders (gen0)
*************************************************************
use "$tmp/final_playonly.dta", clear
foreach d in 4 32 {
	
	qui reghdfe cdrdc`d'_gen0_eth00_a2564 ddtquartile1 $rhs $main
	lincom ddtquartile1
	loc beta1: di %4.3fc r(estimate)
	loc beta = trim("`beta1'")
	loc se1: di %4.3fc r(se)
	loc se = trim("`se1'")
	loc p = r(p)
	reghdfe cdrdc`d'_gen0_eth00_a2564 alpha_* $rhs  $main 
	eventstudy dc`d' 0 `beta' `se' `p' o -30 10 30
	local panel = cond(`d' == 4,"A","B")
	graph export "$out/fig5`panel'.pdf", replace 

}

